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Geologic, chemical and isotopic evidence indicate that Earth has experienced numerous intervals of 
widespread glaciation throughout its history, with roughly 11% of present day Earth's land surface 
covered in ice. Despite the pervasive nature of glacial ice both today and in Earth's past and the 
potential contribution of these systems to global biogeochemical cycles, the composition and 
phylogenetic structure of an active microbial community in subglacial systems has yet to be 
described. Here, using RNA-based approaches, we demonstrate the presence of active and 
endogenous archaeal, bacterial and eukaryal assemblages in cold (0-1 °C) subglacial sediments 
sampled from Robertson Glacier, Alberta, Canada. Patterns in the phylogenetic structure and 
composition of subglacial sediment small subunit (SSU) ribosomal RNA (rRNA) assemblages 
indicate greater diversity and evenness than in glacial surface environments, possibly due to 
facilitative or competitive interactions among populations in the subglacial environment. The 
combination of phylogenetically more even and more diverse assemblages in the subglacial 
environment suggests minimal niche overlap and optimization to capture a wider spectrum of the 
limited nutrients and chemical energy made available from weathering of bedrock minerals. The 
prevalence of SSU rRNA affiliated with lithoautotrophic bacteria, autotrophic methane producing 
archaea and heterotrophic eukarya in the subglacial environment is consistent with this hypothesis 
and suggests an active contribution to the global carbon cycle. Collectively, our findings 
demonstrate that subglacial environments harbor endogenous active ecosystems that have the 
potential to impact global biogeochemical cycles over extended periods of time. 
The ISME Journal (20^3) 7, 1402-1412; doi:10.1038/ismej.2013.31 ; published online 14 March 2013 
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Introduction 

Glacial ice covers ^11% of Earth's land mass today 
and a much larger fraction at times in Earth's past 
(for example, in the late Neoproterozoic, so called 
'Snowball Earth' - 650 Mya; Kirschvink et ah, 2000). 
The advance and retreat of ice during glacial- 
interglacial cycles is thought to have driven large 
scale shifts in the distribution of taxa to habitats 
with more stable climatic conditions (refugia; 
Hodson et ah, 2008) enabling biodiversity to persist. 
Recently, it was suggested that glacial beds, which 
tend to be comparably more stable than other 
surface environments (Skidmore et ah, 2005), may 
have served as refugia for biodiversity during 
periods of inclement climatic, atmospheric or 
geologic conditions (Hodson et ah, 2008). Here, the 
continual exposure of fresh mineral surfaces at the 
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bed of the glacier would serve as a steady source of 
nutrients and energy capable of supporting a 
functionally diverse assemblage of microorganisms 
over extended periods of time (Sharp et ah, 1999; 
Skidmore et ah, 2000, 2005). One-dimensional 
models that integrate metabolic potentials of puta- 
tive microbial assemblages in sub-ice sheet sedi- 
ments during glacial periods indicate that such 
ecosystems could contribute significantly to global 
biogeochemical cycles, specifically the carbon cycle 
(Wadham et ah, 2008, 2012). However, our under- 
standing of life in subglacial sediment environments 
is incomplete and evidence for active microbiomes 
in such systems has yet to be firmly demonstrated. 

Molecular phylogenetic analysis of RNA affords 
the unique opportunity to examine both the compo- 
sition and structure of the active fraction of microbial 
assemblages sampled from natural environments 
(Felske et ah, 1998). Single-stranded RNA is synthe- 
sized only by active cells and it degrades relatively 
rapidly once produced when compared to DNA 
(Hirsch et ah, 2010). Thus, RNA-based approaches 
have been employed to identify the functioning 
members of natural microbial communities. Despite 
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the clear advantage of RNA-based approaches to 
microbial ecology research, such approaches have 
received far less attention in subsurface research. 
This is due primarily to the low biomass associated 
with many subsurface environments (Whitman 
et aL, 1998) including subglacial systems 
(Kastovska et aL, 2007), a feature that is often 
confounded by the presence of complex organic 
and mineral matrices in these systems (Hodson et aL, 
2008) that interfere with the successful extraction 
and purification of intact RNA. However, recent 
technological advances now make it possible to 
characterize assemblages from biomass-limited eco- 
systems using RNA-based approaches, including the 
deep subsurface (Ogram et aL, 1995). These 
advances facilitate not only examination of the 
composition of communities, but also their phylo- 
genetic structure. Patterns in the phylogenetic 
structure of assemblages can provide unique insight 
into the relative roles of stochastic (species neutral) 
and deterministic (competitive or facilitative inter- 
species interactions or environmental filtering) pro- 
cesses in dictating the assembly of macrobiological 
communities (Webb et aL, 2002; Cavender-Bares 
et aL, 2004; Horner-Devine and Bohannan, 2006; 
Bryant et aL, 2008; Cavender-Bares et aL, 2009; 
Meuser et aL, 2013). Such approaches hold tremen- 
dous promise in the study of microbial ecosystems 
where they can provide critical insight into the 
mechanisms that generate and maintain diversity, 
such as speciation, extinction, dispersal and species 
interactions (Martiny et aL, 2006). Nonetheless, they 
have only recently been applied to microbial 
assemblages (Horner-Devine and Bohannan, 2006; 
Bryant et aL, 2008; Meuser et aL, 2013). 

Robertson Glacier (RG), Alberta, Canada is an 
alpine glacier in the Canadian Rockies that, for the 
past decade, has served as a model system to 
understand the contribution of glacial systems to 
local and global biogeochemical cycles (Sharp et aL, 
2002; Boyd et aL, 2010, 2011). Several indirect lines 
of evidence suggest that the subsurface environment 
beneath RG may host an active microbiome. First, 
geochemical data indicate the production of sig- 
nificant amounts of sulfate in the cold (0-1 °C) 
subglacial environment, likely through the weath- 
ering of pyrite present in bedrock (Sharp et aL, 
2002). Second, genetic and microcosm studies 
reveal that RG subglacial sediments host a diversity 
of organisms that are capable of performing several 
redox transformations ex situ, including methano- 
genesis (Boyd et aL, 2010) and nitrate reduction and 
nitrification (Boyd et aL, 2011). Third, the detection 
of coenzyme M, a biomarker of methanogens and 
archaeal methanotrophs (Elias et aL, 1999; Hallam 
et aL, 2004), was detected in subglacial sediments at 
a concentration corresponding to ^3000 cells gram 
dry weight sediment"^ (Boyd et aL, 2010). However, 
it is not known if the populations were active or in a 
state of quiescence because of an incomplete under- 
standing of the stability of coenzyme M under these 



environmental conditions. In this study, we evalu- 
ated the potential for RG subglacial sediments to 
host active microbial communities by quantifying 
and sequencing archaeal, bacterial, and eukaryal 
small subunit (SSU) ribosomal RNA (rRNA) genes 
and their transcripts in genomic DNA and RNA 
extracts, respectively. In order to test for endogeni- 
city in the subglacial ecosystem, these data were 
then compared with data from representative 
surficial glacial environments (snow debris and 
cryoconite sediments) that when melted or dis- 
turbed, could introduce exogenous populations to 
the subglacial system via moulins and crevasses. 
Patterns in the phylogenetic structure and composi- 
tion identified in the respective communities are 
discussed in terms of the role of subglacial systems 
in (i) global biogeochemical cycles and (ii) support- 
ing biodiversity during extended periods of incle- 
ment climatic, atmospheric or geological conditions 
(Hodson et aL, 2008). 



Materials and methods 

Sample collection 

Samples for nucleic acid extraction were collected 
on 14 October 2010 from RG (115°20'W, 50°44'N) in 
Peter Lougheed Provincial Park, Kananaskis County, 
Alberta, Canada. Fine-grained basal sediments were 
collected aseptically from within an ice cave that 
formed at the terminus of the glacier due to 
subglacial discharge at 1200 hours. Subsamples of 
cryoconite sediment and snow, the latter of which 
contained a noticeable amount of debris, were 
collected from the surface of the glacier at 1300 
hours. A detailed description of sample sites and 
sample collection strategy is provided in the 
supporting online material (SOM), SOM Figure 1 
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Figure 1 qRT-PCR of SSU rRNA transcripts (rRNA) and genes 
(rDNA) from subglacial sediment, cryoconite sediment and snow 
debris. Results are presented as the mean of triplicate qPCR 
assays; error bars represent the s.d. of replicates. Copy number is 
normalized to per gram dry weight sediment. Cryo, cryoconite 
sediment; Sub, subglacial sediment. 
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Table 1 Sequencing and diversity statistics for assemblages sampled at Robertson Glacier. SSU cDNA and SSU rDNA assemblages for 
each sample site are separated by a slash 



Diversity statistics 



Mean phylogenetic diversity^ 



Subglacial 



Cryoconite 



Snow 



Predicted coverage of diversity^ 



Subglacial 



Cryoconite 



Snow 



Archaea 
Bacteria 
Eukarya 



0.165/0.151 
0.517/0.527 
0.730/0.599 



0.164/0.068 
0.509/0.520 
0.549/0.569 



0.062/0.024 

0.434/0.450 
0.429/0.482 



0.985/0.985 
0.962/0.937 
0.942/0.955 



0.956/0.993 
0.974/0.871 
0.905/0.788 



0.985/0.985 
0.992/0.943 
0.970/0.900 



Evenness statistics 



Subglacial Cryoconite Snow Subglacial Cryoconite Snow 



Archaea 0.732/1.149 -0.061/0.811 1.713/0.976 0.333/0.760 0.877/0.413 0.913/1.034 

Bacteria - 0.823/ - 1.015 - 0.714/ - 0.845 - 0.369/ - 0.158 - 1.574/ - 2.320 - 1.067/ - 3.269 0.494/ -0.590 

Eukarya - 1.426/ - 0.231 0.408/0.440 0.885/1.060 - 1.697/ - 1.146 -0.234/0.033 0.949/1.546 



Abbreviations: cDNA, complementary DNA; MPD, mean phylogenetic distance; NRI, net relatedness index; NTl, nearest taxon index; 
OTU, operational taxonomic unit; rRNA, ribosomal RNA; SSU, small subunit. 

''MPD, calculated using representative OTUs defined at 97.0% identities. Higher MPD indices are indicative of a greater phylogenetic diversity. 

Boldface values are statistically significant from randomized samples (null model 3) at P<0.10. 

^The percentage of the predicted diversity (defined at 97.0% identities) that was sampled from each assemblage. 

''The NRI and NTl as an indicator of tree-wide and branch tip phylogenetic evenness, respectively, in each assemblage. Increasingly negative NRI 
and NTl are indicative of phylogenetic overdispersion/evenness whereas increasingly positive values are indicative of phylogenetic clustering. 
Boldface values are statistically significant from randomized samples (null model 3) at P<0.10. 



and SOM Table 1. Samples for DNA-based analyses 
(^Ig) were collected in sterile 1.5-ml microcentri- 
fuge tubes with flame-sterilized spatulas and imme- 
diately flash-frozen in a dry ice-ethanol slurry. 
Sediment aliquots ( ^ 1 g) for RNA were collected 
in sterile 2-ml tubes containing 0.5 ml RNALater 
(Qiagen, Valencia, CA, USA) and flash-frozen in a 
dry ice-ethanol slurry. Samples were stored on dry 
ice during transport to the field station and during 
transport back to Montana State University where 
they were stored at - 80 °C until further processed. 
With the exception of snow samples, all samples for 
RNA extraction were stored in RNALater (Qiagen) 
prior to RNA extraction. Snow samples were kept at 
— 80 °C until they were melted (4 °C) and concen- 
trated via centrifugation (14 000xg, 5min., 4 °C) 
immediately prior to molecular-based analyses 
(additional details of snow processing is provided 
in SOM). 



Nucleic acid extraction, quantification, and generation 
of complementary DNA (cDNA) 

DNA extraction and purification was carried with a 
Fast DNA Spin Kit for Soil (MP Biomedicals, Solon, 
OH, USA) as modified from the manufacturer's 
instructions as previously described (Boyd et ah, 
2007b). DNA was extracted in triplicate from three 
independent ^250mg subsamples of wet weight 
sediment in the case of subglacial and cryoconite 
samples or^250mg of debris concentrated from 
melted snow. Equal volumes of each triplicate 
extraction were pooled for further analyses. RNA 
extraction and initial purification was carried 
out using a FastRNA Spin kit for Soil (MP 
Biomedicals). RNA was extracted in triplicate from 



approximately^ 400 mg of wet weight sediment 
(subglacial and cryoconite samples) or^400mg 
of wet weight debris (snow sample). Following 
initial purification, RNA was subjected to DNAse I 
digestion (Roche, Indianapolis, IN, USA) for Ih at 
room temperature (^22 °C). Following digestion, 
RNA was further purified using a High Pure RNA 
Isolation Kit (Roche) and was stored at — 80 °C in a 
solution of 100% ethanol and 0.3 m sodium 
acetate until further processed. The dry solid 
content of the subglacial sediments, cryoconite 
sediments and snow debris used in molecular 
analyses was determined by drying the residual 
material remaining after nucleic acid extraction at 
80 °C for 24 h. 

The concentration of DNA or RNA in each pooled 
volume was determined using a Qubit dsDNA HS 
Assay kit or a Qubit RNA Assay kit (Molecular 
Probes, Eugene, OR, USA) and a Qubit 2.0 Fluo- 
rometer (Invitrogen, Carlsbad, CA, USA), respec- 
tively. As an additional precaution, RNA extracts 
were screened for the presence of contaminating 
genomic DNA by performing a PCR using ^ 1 ng of 
RNA as template and archaeal and bacterial 16S 
rRNA gene primers, described below. Equal volumes 
of each RNA extract were pooled and subjected to 
cDNA synthesis. cDNA was synthesized from 20 ng 
of purified RNA extracted from the subglacial and 
cryoconite sediments and the snow debris using the 
iScript cDNA Synthesis Kit (Bio-Rad Laboratories, 
Hercules, CA, USA) and the following reaction 
cycling conditions: 5min at 25 °C, 30min at 42 °C 
and 5 min at 85 °C. Following synthesis of cDNA, 
samples were purified by ethanol precipitation 
and re-suspended in nuclease-free water for use in 
pyrotag sequencing. 
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PGR amplification of SSU rRNA genes from genomic 
and complementary DNA 

Purified genomic DNA extracts or cDNA were 
subjected to amplification of archaeal and bacterial 
16S ribosomal DNA (rDNA) using primers 344F 
(5'-ACGGGGYGCAGCAGGCGCGA-30 and 915R 
(5'-GTGCTCCCCCGCCAATTCCT-30 at an annealing 
temperature of 61 °C or primers HOOF (5'-YAACGA 
GCGCAACCC-30 and 1492R (5'-GGTTACCTTGTTA 
CGACTT-30 and an annealing temperature of 55 °C, 
respectively, (Boyd et al., 2007a). For amplification 
of eukaryotic 18S rDNA, primers euk-A7F (5'-AACC 
TGGTTGATCCTGCCAGT-30 (Medlin et al, 1988) 
and Euk-570R (5'-GCTATTGGAGCTGGAATTAC-30 
(Weekers et al., 1994) were used at an annealing 
temperature of 42 °C. For each set of primers, ^ 1 ng 
of purified genomic DNA or cDNA was subjected to 
PGR in triplicate using the following reaction 
conditions: initial denaturation at 94 °C (4min), 
followed by 35 cycles of denaturation at 94 °C 
(Imin), annealing at the optimal temperature for 
each primer pair (Imin), primer extension at 72 °C 
(1.5 min), followed by a final extension step at 72 °C 
for 20 min. Reactions contained 2mM MgCls (Invi- 
trogen), 200 jiM each deoxynucleotide triphosphate 
(Eppendorf, Hamburg, Germany), 0.5 jiM forward 
and reverse primer (Integrated DNA Technologies, 
Coralville, lA, USA), 0.4 mg ml"^ molecular-grade 
bovine serum albumin (Roche), and 0.25 units Taq 
DNA polymerase (Invitrogen) in a final reaction 
volume of 50|il. Positive control reactions were 
performed using genomic DNA from Azotobacter 
vinelandii DJ, Roseiflexus castenholzii, Sulfolobus 
solfataricus P2 and Chlamydomonas reinhardtii. 
Negative control reactions were performed in the 
absence of added genomic DNA or cDNA. 



Quantitative PGR (qPGR) 

qPCR was used to estimate the number of SSU rDNA 
templates in genomic DNA extracted from subglacial 
sediments, cryoconite sediment and snow debris. 
qPCR followed protocols described previously 
(Boyd et al, 2011). Briefly, SSU rDNA genes 
amplified from the archaeal, bacterial and eukaryal 
positive controls (as listed in previous section) were 
used to generate standard curves for use in relating 
template copy number to threshold qPCR amplifica- 
tion signal. Triplicate PGR products were purified 
using a QIAquick PGR Purification Kit (Qiagen) and 
quantified as described above. The concentrations of 
each replicate PGR product varied by less than a 
factor of 1.3 and thus were averaged for use in 
generating standard curves. An archaeal 16S rDNA 
standard curve was generated over 5 orders of 
magnitude from 2.6 x 10^ to 2.2 x 10^ copies of 
template per assay (R^ = 0.998). A bacterial 16S 
rDNA standard curve was generated over 5 orders of 
magnitude from 3.1 x 10^ to 5.7 x 10^ copies of 
template per assay (R^ = 0.998). An eukaryotic 18S 
rDNA standard curve was generated over 5 orders of 
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magnitude from 9.1 x 10^ to 3.7 x 10^ copies of 
template per assay (R^ = 0.997). The detection limit 
for all assays was between 10 and 15 copies of 
template per assay. qPGR assays were performed in 
optically clear PGR tubes (Qiagen) in a Rotor-Gene 
300 quantitative real-time PGR machine (Qiagen) 
using a SsoFast EvaGreen Supermix qPGR Kit (Bio- 
Rad Laboratories). Assay reactions were amended 
with molecular-grade bovine serum albumin to a 
final concentration of 0.4 mg ml"^ (Roche). qPGR 
cycling conditions were as follows: initial denatura- 
tion (95 °G for 10 min) followed by 40 cycles of 
denaturation (95 °G for 10 s), annealing (at the 
optimal annealing temperature for each primer as 
described above for 10 s) and extension (72 °G for 
20 s). Specificity of the qPGR assays was verified by 
melt curve analysis. The reported template abun- 
dances are the average and standard deviation of 
qPGR assays performed in triplicate. 



Quantitative reverse transcription-PGR (qRT-PGR) 
qRT-PGR was used to quantify the abundance of 16S 
rRNA transcripts in RNA extracted from subglacial 
sediment, cryoconite sediment and snow debris. 
qRT-PGR was performed with the Power SYBR 
Green RNA-to-Gx one-step kit from Invitrogen 
according to the manufacturer's protocol and were 
assayed on a RotorGene-Q real-time PGR detection 
system from Qiagen. Reactions were performed in 
triplicate with each reaction containing ^lOng of 
total RNA, 500 UM forward and reverse primer, in a 
final reaction volume of 20 jil. The following cycling 
conditions were employed: reverse transcription at 
48 °G (30 min) followed by initial activation of the 
DNA polymerase at 95 °G (10 min) followed by 40 
cycles of denaturation at 95 °G (15 s), annealing and 
extension (at the optimal annealing temperature for 
each primer as described above for Imin). Specifi- 
city of the qRT-PGR assays was verified by melt 
curve analysis. Gontrol reactions to ensure that the 
amplified product was not due to genomic DNA 
contamination were performed in the absence of 
reverse transcriptase or in the absence of added 
template RNA. An archaeal 16S rDNA standard 
curve was generated over 5 orders of magnitude 
from 6.5 x 10^ to 3.5 x 10^ copies of template per 
assay (R^ = 0.997). A bacterial 16S rDNA standard 
curve was generated over 5 orders of magnitude 
from 9.7 x 10^ to 2.3 x 10^ copies of template per 
assay (R^ = 0.998). A eukaryotic 18S rDNA standard 
curve was generated over 5 orders of magnitude 
from 4.4 x 10^ to 5.1 x 10^ copies of template per 
assay (R^ = 0.997). The detection limit for all assays 
was between 14 and 18 copies of template per assay. 



Sequence analysis 

SSU rDNA and SSU cDNA were sequenced by the 
Research and Testing Laboratory (Lubbock, TX, 
USA). SSU rDNA and cDNA from each fraction 
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were barcoded as previously described (Dowd et ah, 

2008) using the primers described above for SSU 
rDNA analysis and were sequenced using a 454 
Genome Sequencer FLX System (Roche, Nutley, 
NJ, USA). Each sample was sequenced once. Post 
sequence processing was performed using the Mothur 
(ver. 1.24.1) sequence analysis platform (Schloss et ah, 

2009) . Raw libraries were trimmed to a minimum 
length of 200 bases and were subjected to a filtering 
step using the quality scores file to remove sequences 
with anomalous base calls. Unique sequences were 
aligned using SILVA databases specific for SSU rDNA 
from each domain and sequences were trimmed using 
a defined start and end site based on inclusion of 75% 
of the total sequences; sequences that started before or 
after these defined positions were removed without 
further consideration. The resulting unique sequences 
were pre-clustered to remove amplification and 
sequencing errors and chimeras were identified and 
removed using UCHIME (Edgar et ah, 2011). Opera- 
tional taxonomic units (OTUs) were assigned at a 
sequence similarity of 97.0% using the furthest- 
neighbor method. The remaining sequences were 
randomly sub-sampled in order to normalize the total 
number of sequences in each library. Collectively, 
these steps resulted in a normalized size of 1136, 1278 
and 938 sequences for each archaeal, bacterial and 
eukaryal library, respectively. Rarefaction curves were 
used to compute the percent coverage of the predicted 
taxonomic richness for each library, which indicating 
that >87% of the predicted sequence richness was 
sampled for each library (Table 1). Sequences were 
classified using the Bayesian classifier (Wang et ah, 
2007) and the RDP database, with manual verification 
using BLASTn (Supplementary Tables 2-A). Sequences 
representing each OTU have been deposited in the 
NCBI SRA database under accession numbers 
SAMNOl 729343 to SAMN01730044 and KC5 72493 
(Supplementary Tables 2-A). 

Sequences representing each unique OTU (defined 
at 97.0% sequence identities) were compiled for 
each domain. ClustalX (ver. 2.0.9; Larkin et aL, 2007) 
was used to align nucleic acid sequences using 
default parameters. Archaeal, bacterial and eukaryal 
SSU rDNA and SSU cDNA alignments were 
subjected to evolutionary model prediction via 
jModeltest (ver. 2.1.1; Darriba et aL, 2012), Max- 
imum-likelihood phylogenetic reconstruction via 
PhyML (version 3.0; Guindon and Gascuel, 2003) 
specifying the general time reversible model and 
gamma distributed rate variation with a proportion 
of invariable sites, and rate smoothing using the 
multidimensional version of Rambaut's parameter- 
ization as implemented in PAUP (ver. 4.0; Swofford, 
2001) as previously described (Meuser et aL, 2013). 



Community ecology analysis 

Rate smoothed ultrameric trees were used to 
calculate the mean phylogenetic distance (MPD), 
net relatedness index (NRI) and nearest taxon index 



(NTI) for each assemblage while specifying among 
individual abundances (-a parameter) over 999 
iterations with the program Phylocom (ver. 4.0.1; 
Webb et aL, 2008). MPD is an abundance weighted 
metric that describes the pairwise phylogenetic 
distance between sequences in a community, when 
compared with the total sequence pool. Assem- 
blages with higher MPD indices exhibit a greater 
phylogenetic diversity relative to the total sequence 
pool. The NRI is a measure of tree-wide phyloge- 
netic clustering of sequences, whereas NTI is more 
sensitive to clustering at the terminals of the tree 
(Webb et aL, 2002). Increasingly positive NRI and 
NTI metrics indicate that co-occurring species are 
more phylogenetically related than expected by 
chance (phylogenetic clustering). In contrast, 
increasingly negative NRI and NTI metrics indicate 
that co-occurring species are less phylogenetically 
related than expected by chance (phylogenetic 
overdispersion). We tested whether these values 
differed significantly from that of a randomly 
assembled community by comparison with a null 
model generated by the independent swap method 
(null model 3 within Phylocom). The independent 
swap method was selected because of the results of a 
previous empirical study (Kembel, 2009) that 
demonstrated that of the models available in 
Phylocom, this method is the most robust for 
detecting niche-based assembly processes and is 
less prone to type I errors. One thousand swaps were 
performed over 999 permutations, specifying 
sequence abundance weights (-a) in the analysis. A 
two-tailed significance test was used to evaluate the 
rank of observed values. When >900 permutations 
supported the observed values rather than the 
random or null model (P<0.10), the observed rank 
was deemed to be significant. 

Phylocom was also used to calculate Rao's com- 
munity phylogenetic relatedness for archaeal, bac- 
terial and eukaryal assemblages using relative 
sequence abundance weights and rate-smoothed 
ultrameric trees. PAST (ver. 2.17b; Hammer et aL, 
2001) was used to generate cluster dendograms 
specifying the single linkage method and Euclidean 
distances. Bootstrap values correspond to the fre- 
quency that each node was observed in a given 
position out of 1000 replicates. 



Results and discussion 

Abundance of SSU rRNA and rDNA 
The detection of archaeal, bacterial and eukaryal 
SSU rRNA in subglacial sediments suitable for 
reverse transcription to cDNA (Figure 1) provided 
the first direct evidence that the RG subglacial 
system harbors active microbial populations. Sur- 
prisingly, subglacial sediments harbored only a 
slightly lower total (archaea + bacteria + eukarya) 
SSU rDNA gene abundance than surface snow or 
cryoconite (Figure 1), although the total SSU rRNA 
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abundance was lower in the subglacial sediments 
than surface snow or cryoconite. As such, the ratio 
of total SSU rRNA to SSU rDNA genes was 
significantly lower (P<0.03 in all two-tailed pair- 
wise T-tests) in the subglacial sediments 
(0.22 + 0.01) than in cryoconite debris and surface 
snow (0.42 + 0.03 and 0.27 + 0.02, respectively; 
Figure 1). The lower abundance of SSU rRNA and 
lower ratio of SSU rRNA to SSU rDNA in the 
subglacial sediments could be attributable to the 
dependence on limited chemical energy and oligo- 
trophic conditions in the RG subglacial environment 
(Boyd et ah, 2011) when compared with the surface 
environment, which has direct access to light energy 
capable of supporting photosynthesis. 



Phylogenetic structure and diversity of assemblages 
Unlike traditional alpha diversity metrics (for 
example, pairwise comparison of OTU overlap 
among communities), beta diversity metrics such 
as NRI, NTI and MPD use a phylogenetic framework 
for determining the extent of evolutionary diver- 
gence between taxa that comprise a community, 
which can vary considerably (Lozupone and Knight, 
2008). Given that many metabolic traits are con- 
served during the evolution of a lineage (Blomberg 
et aL, 2003), positive relationships between the 
phylogenetic relatedness of species and their overall 
ecological similarity are often observed (Webb et aL, 
2002; Wiens and Graham, 2005; Cadotte et aL, 2009). 
For this reason, beta metrics of diversity (for 
example, MPD) can provide insight into both the 
taxonomic and ecological variation among commu- 
nities. Similarly, beta metrics of phylogenetic struc- 
ture, which describe the extent to which phylotypes 
comprising a community are clustered on ultrameric 
phylogenetic trees (for example, NRI/NTI), provide 
insight into the role of stochastic (species neutral) 
and deterministic mechanisms, either biological 
(competitive interactions) or environmental (habitat 
filtering), in the assembly of the community (Webb 
et aL, 2002; Cavender-Bares et aL, 2004; Horner- 
Devine and Bohannan, 2006). This information can 
provide a basis for predicting a number of properties 
of an ecosystem, including overall productivity 
given finite resources (Maherali and Klironomos, 
2007), the extent of niche overlap (Cavender-Bares 
et aL, 2004), as well as the susceptibility of a 
community when confronted with environmental 
change (Knapp et aL, 2008). 

Archaeal assemblages at RG tended to be phylo- 
genetically clustered as evinced by positive NRI 
and/or NTI metrics regardless of whether they were 
derived from genomic DNA or cDNA fractions or if 
they were derived from surface cryoconite, surface 
snow or subglacial sediment. This suggests that 
archaea dispersing into each of the respective 
environments are excluded or included based on 
their ecologies, which is dictated largely by the 
physiological traits harbored by an organism (Webb 
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et aL, 2002; Horner-Devine and Bohannan, 2006; 
Martiny et aL, 2006). In other words, the conditions 
of the environment, which are the combined result 
of abiotic and biotic factors, are likely to select for a 
particular set of traits that are compatible with the 
characteristics of an environment and which permit 
colonization. The inferred physiologies of archaeal 
populations (described in detail below) in the 
subglacial and cryoconite assemblages is that of 
strict anaerobes, whereas those in the surface snow 
assemblages are strict aerobes. Thus, aerobic/anae- 
robic boundaries may select for specific traits related 
to the utilization or avoidance of oxygen in the 
respective ecosystems. Such a situation is consistent 
with data that suggest anoxic conditions prevail at 
the bed of an alpine glacier in Switzerland (Tranter 
et aL, 2002). 

In contrast, bacterial and eukaryal SSU cDNA 
assemblages associated with subglacial sediments 
and cryoconite exhibited significant and negative 
NRI and/or NTI, with those associated with the 
subglacial sediments typically exhibiting more 
statistically significant and negative values 
(Table 1). Increasingly negative NRI/NTI indices 
suggest that assemblages contain phylotypes that are 
more phylogenetically divergent or that are affiliated 
with higher taxonomic ranks, given the species 
richness of the system (that is, phylogenetically 
overdispersed or phylogenetically even) (Webb 
et aL, 2002). Evidence for evenly structured com- 
munities has been interpreted to be the result of (i) 
competitive interactions among closely related 
species that share a fundamental niche thereby 
limiting their co-existence over the long term 
(competitive exclusion) (Losos 1994; Webb et aL, 
2002; Cavender-Bares et aL, 2004; MaheraU and 
Klironomos, 2007) or (ii) facilitative interactions 
where facilitator species create microhabitats that 
support distantly related specie(s) with non-over- 
lapping niches (Lortie, 2007; Valiente-Banuet and 
Verdu, 2007). Evenly structured soil communities 
are associated with the utilization of a greater 
spectrum of the resources when compared with 
phylogenetically clustered assemblages (Maherali 
and Klironomos, 2007; Fornara and Tilman, 2008). 
Thus, phylogenetically even assemblages in the RG 
subglacial environment may reflect the oligotrophic 
conditions present in the RG subsurface environ- 
ment (Boyd et aL, 2010, 2011) and selection to 
optimize the utilization of limited resources. 

In this study, the MPD of archaeal, bacterial and 
eukaryal SSU cDNA was higher in the subglacial 
ecosystem than surface snow or cryoconite debris 
ecosystems (Table 1), consistent with conclusions 
drawn from analyses of NRI/NTI metrics. As closely 
related species tend to have similar physiological 
traits and often have overlapping ecological niches 
(Webb et aL, 2002; Wiens and Graham, 2005), it 
follows that distantly related species should have 
divergent physiological traits, a feature that would 
help maintain minimal niche overlap and 
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potentially minimize competitive interactions 
(Maherali and Klironomos, 2007). Alternatively, 
the tendency for overdispersed subglacial assem- 
blages to exhibit elevated phylogenetic diversity 
may point to an important role for facilitation among 
divergent, co-occurring species. Previous studies 
have shown that facilitative or mutualistic inter- 
actions among plant populations increase the 
availability of resources and alter substrate char- 
acteristics (Callaway, 1995), features that may lead 
to increased phylogenetic diversity (Valiente-Banuet 
and Verdu, 2007) and which may stabilize ecosys- 
tems (Ringel et aL, 1996). Mutualistic relationships 
are common in nature and may be more pronounced 
in nutrient-limited ecosystems such as RG, which 
has been shown to be oligotrophic with respect to 
bioavailable nitrogen (Boyd et aL, 2011) and phos- 
phorus (ES Boyd, unpublished). 

The greater phylogenetic diversity of active (RNA- 
based SSU cDNA analysis) and total (DNA-based 
SSU rDNA analysis) archaea, bacteria and eukarya 
associated with the subglacial sediments when 
compared with surface assemblages was unantici- 
pated. In particular, the detection of active eukarya 
in the hypoxic subglacial sediments was surprising 
given the lack of previous evidence for this domain 
in cold subglacial environments. Recent evidence, 
however, suggests that eukaryote populations are 
active in a number of deep sea marine sediment 
environments characterized by cold, hypoxic and 
oligotrophic conditions (Edgcomb et aL, 2010; 
Schippers et aL, 2012). In deep marine sediments, 
eukaryotes are likely to be bacterivores that graze 
bacterial cells or heterotrophs dependent on buried 
organics. Phylogenetic inference of the dominant 
active eukaryotes in RG subglacial sediments (see 
below) indicate that they are heterotrophic and 
dependent on organic carbon. This finding is 
consistent with the abundance of particulate organic 
carbon (25±14mg g sediment"^; Boyd et aL, 2010) 
and pore water dissolved organic carbon (60 jimol 1" ^; 
Boyd et aL, 2011) in RG sediments. 



Comparison of SSU rDNA and cDNA composition 
Archaeal, bacterial and eukaryal SSU cDNA and 
rDNA assemblages from each environment tended to 
form sister groups or branch closely together in 
cluster analysis (Figure 2), providing evidence that 
the SSU cDNA in each subglacial assemblage is the 
result of endogenous synthesis. In the case of each 
domain, the subglacial and cryoconite assemblages 
clustered closely together and distantly from the 
snow assemblage, suggesting that the subglacial 
sediment and cryoconite ecosystems are more 
similar ecologically. Intriguingly, the depth (or the 
extent to which the assemblages differ) of the 
clusters that comprise the SSU cDNA and rDNA 
sister groups for the subglacial and cryoconite 
archaeal, bacterial and eukaryal microbiomes was 
generally shallower than that from the snow 
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Figure 2 Cluster dendograms depicting the Rao phylogenetic 
relatedness of archaeal (a), bacterial (b) and eukaryal (c) SSU 
rDNA and SSU cDNA assemblages. The cophenetic correlation 
coefficient indicating the extent to which the pairwise distances 
between the original unmodeled dissimilarity data points 
were preserved in each dendogram was 0.9997, 0.9486 and 
0.8791, respectively. Cryo, cryoconite sediment; Sub, subglacial 
sediment. 



microbiome, suggesting that the majority of the 
SSU rDNA genes present in the subglacial and 
cryoconite microbiomes are being transcribed, with 
a much lower fraction of the SSU rDNA genes being 
transcribed in the snow microbiome at the time of 
sampling (midday). 



Composition of archaeal, bacterial and eukaryal 
assemblages 

The composition of SSU rDNA and cDNA commu- 
nities was similar in each respective environment 
especially in the subglacial ecosystem (Figure 3), 
consistent with the results of the cluster analysis. 
In general, the inferred physiology of the dominant 
populations associated with subglacial sediments 
suggests the presence of a complex ecosystem that 
contains both aerobic and anaerobic taxa 
likely harbored in numerous microhabitats 
(Supplementary Tables 2-4). Archaeal SSU cDNA 
was dominated (70% of total) by sequences 
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Figure 3 Composition of SSU rDNA and SSU cDNA in nucleic 
acids recovered from subglacial sediment, cryoconite sediment 
and snow debris. Representative OTUs for each library were 
binned at the order level for eukarya (a) and bacteria (b) and at the 
family level for archaea (c). Taxonomic bins that represen- 
ted<5.0% of the total sequences from each assemblage were 
pooled and depicted as 'Other'. Cryo, cryoconite sediment; Sub, 
subglacial sediment. 



affiliated with the strictly anaerobic methanogen 
order Methanomicrobiales with a smaller fraction 
(25% of total) affiliated with the Methanosarcinales 
(Supplementary Table 2). In contrast, bacterial SSU 
cDNA was dominated by sequences affiliated with 
aerobic autotrophs Sideroxydans sp. (12% of total), 
Sulfurihydrogenibium sp. (12% of total) and Piano- 
coccus sp. (10% of total) (Supplementary Table 3). 
Sequences affiliated with Sulfurihydrogenibium sp., 
which were only identified in cold (0-1 °C) sub- 
glacial SSU cDNA and rDNA assemblages, are of 
particular interest as sequences affiliated with this 
lineage have yet to be identified in any environment 
with a temperature of <50°C (Reysenbach et al., 
2005). The eukaryal SSU cDNA assemblage in the 
subglacial ecosystem was dominated by sequences 
affiliated with two unicellular orders, including the 
ciliate order Strichotrichida [Amphisiella sp., 29% of 
total) and the amoeba order Tectofilosida [Capsellina 
sp., 17% of sequences), both of which are hetero- 
trophic (Supplementary Table 4). The eukaryal 
assemblage associated with the cDNA fraction 
exhibited differences from that associated with the 
DNA fraction, including a greater abundance of 
sequences affiliated with the amoeba Stichotrichida 
and the flagellate Thaumatomonadida and, impor- 
tantly, an absence of sequences affiliated with the 
green alga Ulotrichales (Figure 3; Supplementary 
Table 4). The recovery of active and dominant 
archaea, bacteria and eukarya in subglacial sedi- 
ments that are putatively involved in autotrophic 
methane production, lithoautotrophic carbon diox- 
ide fixation and heterotrophy, respectively, suggests 
that these ecosystems have significant potential to 
impact global biogeochemical cycles, in particular 
the carbon cycle. These results confirm previous data 
that were suggestive of the presence of an active 
methanogen assemblage in this environment (Boyd 
et al., 2010) and corroborate recent reports suggesting 
that these ecosystems may contribute to the global 
methane cycle over recent glacial-interglacial cycles 
(Wadham et al, 2008, 2012). 

Comparison of the taxonomic affiliation of SSU 
cDNA derived from subglacial sediment, surface 
snow and cryoconite communities reveal distinct 
differences among the communities (Figure 3; 
Supplementary Tables 1-3). In general, the 
cryoconite and subglacial sediment SSU cDNA 
communities were most similar, with pronounced 
differences noted between these assemblages and 
the assemblages associated with snow. In particular, 
the abundance of eukaryal sequences affiliated with 
Chlamydomonas sp., an alga, in the surface snow 
(67% of total) and cryoconite (40% of total) but not 
the subglacial sediment is consistent with previous 
hypotheses suggesting that subglacial assemblages 
are supported by chemical energy (Skidmore et al., 
2000, 2005; Gaidos et al., 2008) while surface 
communities are supported primarily by light 
energy (Telling et al., 2010). Sequences affiliated 
with the green alga Ulotrichales were detected in 
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both the SSU cDNA and rDNA fraction of the 
cryoconite assemblage, but only in the rDNA 
fraction of the subglacial sediment assemblage. This 
may indicate the subglacial environment receives 
exogenous material from the surface, such has been 
demonstrated in a separate Canadian Arctic glacial 
system using DNA-based fingerprinting methods 
(Bhatia et aL, 2006). However, the physical char- 
acteristics of the environment (for example, lack of 
available light) are likely to preclude activity of this 
phototrophic population, which may explain the 
lack of detectable sequences affiliated with this 
taxon in the cDNA fraction. Such observations 
further support endogenous synthesis of RNA in 
the subglacial environment. 

Subglacial sediments and cryoconite harbored 
similar archaeal SSU cDNA assemblages, both of 
which were dominated by sequences affiliated with 
methanogens within the Methanomicrobiales and to a 
lesser extent Methanosarcinales (Figure 3). In con- 
trast, the surface snow archaeal SSU cDNA assem- 
blage was dominated (98% of total) by sequences 
affiliated with nitrifying archaea (Nitrososphaerales). 
Bacterial SSU cDNA assemblages also exhibited 
distinct differences, most notably the detection of 
sequences affiliated with Aquificales only in the 
subglacial community, and the detection of sequences 
affiliated with the iron utilizing Gallionellales 
[Sideroxydans sp.) only in the subglacial sediment 
and cryoconite debris communities. 



Conclusions 

The data presented here provide the first RNA-based 
evidence for an active, diverse and endogenous 
microbial community adapted to the dark and 
oligotrophic conditions that characterize subglacial 
environments. At RG, the subglacial sediments 
appear to host phylogenetically more diverse assem- 
blages than cryoconite and surface snow assem- 
blages. This observation, coupled with phylogenetic 
evidence that is suggestive of competitive or facil- 
itative interactions among members of the subglacial 
microbiome, may indicate that the assembly of the 
community reflects the limited nutrient and energy 
condition of the subglacial system. We propose that 
the continual exposure of fresh mineral surfaces due 
to bedrock comminution by flowing glacial ice 
provides ample mineral-based energy and nutrients 
capable of sustaining an active and stable subglacial 
ecosystem during extended periods of widespread 
glaciation (Snowball Earth, for example; Hoffman 
et al, 1998; Kirschvink et aL, 2000). 
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